#ifndef RECONSTRUCTION_PRECALC_DIFFS_HPP
#define RECONSTRUCTION_PRECALC_DIFFS_HPP

namespace hms::reconstruction
{

template<typename LR_w, typename LR_r, typename LR_diff, typename Dx>
void musclScheme(
	LR_w& phiLR,
	LR_w& psiLR,
	const LR_r& phiCellLR,
	const LR_diff& phiDiffL,
	const LR_diff& phiDiffM,
	const LR_diff& phiDiffR,
	const Dx& dxQuotientL,
	const Dx& dxQuotientR,
	const LimiterFunction& limFunc
){
	using namespace reconstruction;
	constexpr Index nRows { nRowsLR<LR_w>() };

	/* we write rLR to phiLR because the size fits */
	LR_w& rLR {phiLR};

	/* with precalculated diffs */
	for ( int i{0}; i<2; ++i ){
		side(rLR, i) = select(
			phiDiffM == 0, 0,
			(i == 1 ? phiDiffR : phiDiffL ) / phiDiffM
			* repl<nRows>( i == 1 ? dxQuotientR : dxQuotientL )
		);
	}

	/* map to the same memory space, but as one array, to simplify access */
	auto&&
		psi { mapOrRefLR<ArrayNXsMap<nRows>>(psiLR) },
		r   { mapOrRefLR<ArrayNXsMap<nRows>>(rLR  ) };

	setSlopeLimiter( psi, r, limFunc );

	for (int i{0}; i<2; ++i){
		side(phiLR,i) =
			side(phiCellLR,i) + 0.5 * side(psiLR,i) *
			phiDiffM * (i == 1 ? -1 : 1);
	}
}

template<typename SolverType,
	typename EV,
	typename LR_r1, typename LR_r2,
	typename LR_diff1, typename LR_diff2,
	typename Dx>
void highOrder(
	EV& ev,
	const LR_r1& dCellLR,
	const LR_diff1& dDiffL,
	const LR_diff1& dDiffM,
	const LR_diff1& dDiffR,
	const LR_r1& hCellLR,
	const LR_diff1& hDiffL,
	const LR_diff1& hDiffM,
	const LR_diff1& hDiffR,
	const LR_r2& qCellLR,
	const LR_diff2& qDiffL,
	const LR_diff2& qDiffM,
	const LR_diff2& qDiffR,
	const Dx& dxL,
	const Dx& dxR,
	const LimiterFunction& limFunc
){
	using S = SolverType;
	using reconstruction::musclScheme;

	using Map1 = Array1XsMap;
	using Map2 = Array2XsMap;
	using LR_w1 = decltype(ev.dLR);
	using LR_w2 = decltype(ev.vLR);
	MapOrRef<Map1, LR_w1>
		d { mapOrRefLR<Map1>(ev.dLR) },
		z { mapOrRefLR<Map1>(ev.zLR) },
		h { mapOrRefLR<Map1>(ev.hLR) };
	MapOrRef<Map2, LR_w2>
		v { mapOrRefLR<Map2>(ev.vLR) };

	/* we use the same memory space for q as for v, because the size is
	 * always the same and all operations are coeff-wise */
	auto& q {v};
	auto& qLR { ev.vLR };

	musclScheme(ev.dLR,ev.psiDH, dCellLR,dDiffL,dDiffM,dDiffR, dxL,dxR,limFunc);
	musclScheme(ev.hLR,ev.psiDH, hCellLR,hDiffL,hDiffM,hDiffR, dxL,dxR,limFunc);
	musclScheme(   qLR,ev.psiQ , qCellLR,qDiffL,qDiffM,qDiffR, dxL,dxR,limFunc);

	z = h - d;
	v = colwiseSelect( d < S::minWaterDepth, 0, q.rowwise() / d );
	d = select( d < S::minWaterDepth, 0, d );
	h = select( d < S::minWaterDepth, z, h );
}


/* provide chaining overloads for compatibility */

#define PRECALC_MUSCL_ARGS(CELL_LR, CELL_LLRR) \
	CELL_LR, \
	left (CELL_LR  ) - left (CELL_LLRR), \
	right(CELL_LR  ) - left (CELL_LR  ), \
	right(CELL_LLRR) - right(CELL_LR  )

template<typename LR_w, typename LR_r, typename Dx>
void musclScheme(
	LR_w& phiLR,
	LR_w& psiLR,
	const LR_r& phiCellLR,
	const LR_r& phiCellLLRR,
	const Dx& dxQuotientL,
	const Dx& dxQuotientR,
	const LimiterFunction& limFunc
){
	musclScheme(
		phiLR, psiLR,
		PRECALC_MUSCL_ARGS(phiCellLR, phiCellLLRR),
		dxQuotientL, dxQuotientR, limFunc
	);
}

template<typename SolverType,
	typename EV,
	typename LR_r1, typename LR_r2,
	typename Dx>
void highOrder(
	EV& ev,
	const LR_r1& dCellLR,
	const LR_r1& dCellLLRR,
	const LR_r1& hCellLR,
	const LR_r1& hCellLLRR,
	const LR_r2& qCellLR,
	const LR_r2& qCellLLRR,
	const Dx& dxL,
	const Dx& dxR,
	const LimiterFunction& limFunc
){
	highOrder<SolverType>( ev,
		PRECALC_MUSCL_ARGS(dCellLR, dCellLLRR),
		PRECALC_MUSCL_ARGS(hCellLR, hCellLLRR),
		PRECALC_MUSCL_ARGS(qCellLR, qCellLLRR),
		dxL, dxR, limFunc
	);
}

#undef PRECALC_MUSCL_ARGS

} // namespace hms::reconstruction

#endif